An Integrable Model for Density-Modulated Quantum Condensates 
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An integrable model possessing inhomogeneous ground states is proposed. The model is a higher-order ana- 
log of the nonlinear Schrodinger equation. The inverse scattering theory with the elliptic-functional background 
is formulated, and the W-soliton solution is obtained as a particular solution. We provide the exact bosonic 
and fermionic quasiparticle eigenstates and reveal various kinds of dynamics such as dark soliton billiards, 
dislocations, gray solitons, and envelope solitons. 



Introduction. — Spatially inhomogeneous quantum conden- 
sates have been attracting a lot of attention for a long time. 
For bosonic condensates, the supersolid phase, which was 
originally discussed four decades ago ||l|-l3|], has received a 
renewed interest since the torsional oscillator experiments of 
'^He ly, ISfl. While the most recent work ||6[] has concluded 
the absence of supersolidity, the candidate of supersolid is 
now also proposed in Rydberg matters 17|,|8|]. For fermionic 
condensates, the realization and observation of Fulde-Ferrell 
(FF) 1911 and Larkin-Ovchinnikov (LO) jlOll states have been 
a long-standing topic. Within a framework of self-consistent 
Bogoliubov-de Gennes (BdG) formalism, the LO state is 
shown to be a ground state in the presence of a magnetic field 



or a population imbalance 111 
are, for example 



12II . Experimental candidates 

CeCoIns in condensed matters 1131 Il4ll 

and the spin-imbalanced superfluid *Li in ultracold atoms 

iIMl. 

To study the quantum condensates, the nonlinear 
Schrodinger (NLS) type equation and its generalizations 
are often used, and referred to as the Gross-Pitaevskii or the 
Ginzburg-Landau (GL) equation for bosonic or fermionic 
systems. Though Gorkov's original derivation ensures that 
the GL description is good only near T - Tc, the recent 
studies show that the gap function obeys the NLS equation 
with higher-order corrections even near T - 1181420(1 . 

Many theoretical studies have established a common and 
model-independent understanding for the mechanism of 
modulation, the properties of ground states, and low-energy 
excitations around them. Compared to stationary states, 
however, the nature of nonlinear excitations such as solitons 
or vortices passing through modulated condensates has not 
been fully investigated yet, because of the difficulty to treat 
time-dependent phenomena. The knowledge for solitons is 
also important to understand transport phenomena past an 
obstacle in non-stationary regimes li21il . To investigate these 
issues, an integrable model would play a prominent role, 
since we can access the various kinds of dynamics exactly. In 
this Letter we propose an integrable model of non-uniform 
quantum condensates. Solving it by the inverse scattering 
method (ISM), we obtain an A^-soliton solution and reveal its 
behavior. We also provide eigenstates of the BdG operators 



for both bosonic and fermionic systems, which are essential 
for the physics of quasiparticles and Nambu-Goldstone (NG) 
modes. 

Note that the solitons given here are different from gap 
solitons. The system forms a pattern not by a periodic external 
force but by itself, and hence the modulated background and 
the solitons influence each other. 

The energy functional of the model proposed in this Letter 
is given by 



// = C3/3 -hcs/s. 



(1) 



where I3 := jdx(}i//J^ + \ip\'^) and I5 := Jdxjli/r^p + 
6|iA|2|i/r,|2 + [(|(/r|2),]2 + 2|i/r|''} are the third and fifth con- 



served quantities in the NLS hierarchy 11221 12311 . How to 
find this model as a modulated condensate is explained be- 
low. The corresponding partial differential equation is given 
by id,iff = 6H/dtfr*. For convenience, we henceforth consider 
the equation for q :- i/'e''". The equation is then given by 

iq, ^- iiq-\- C3{-q^_, -^■ 2\qf-q) -H cs[q^_^j,,.c 

-2(\q\\,q-3q*(q\., + 6\q\''q], (2) 

where the subscripts t and x denote the differentiation. // is a 
chemical potential if 1^ is regarded as a bosonic condensate. 
This is a special case of the AKNS3 equation ll24ll . 

The idea of model construction. — Before beginning the 
analysis of Eq. (O, let us see how to construct an integrable 
model of modulated condensates. We first show that the 
model with a non-local interaction such as soft-core bosons 



1 25142911 , which are used as a model of supersolid, can be ap- 
proximated by a higher-order differential equation. Let us 
demonstrate it by using the Gaussian-like two-body interac- 
tion V{x) - --£i£ie"-*^'^"^'''\ where Vq^x) is a slowly-varying 



even function. If a 

,2 



0, we can show 



1 __Q-xl{Aa'-) 



2a -\Jn 

S{x) + a^6"(x) + y5""(x) + ■ ■ ■ . Using this, the NLS equa- 
tion for the soft-core boson model \dt^{x, t) - -d^.tfr{x, t) H- 
j dyV(x - y)\i(/(y, t)f-if/(x, t) can be approximated up to 0{a^) 
as follows: 

i«A/ = - (Ax. + [VqW\- + milf\\, + V4m\.x.\ iff, (3) 



where % = Vo{0) + a^Vo (0) + ^Vi;"{Q), % = a^iVoiO) + 
3a^V'^{Qi)), and V4 = ^Vo(O). Even though Eq. Q is a too 
rough approximation for the original non-local model, it ex- 
hibits a roton minimum in the Bogoliubov spectrum and has 
an inhomogeneous ground state under certain parameter re- 
gions, as similar to Ref. [27]. So it can be used as a model 
of supersolid. It is reasonable that the higher-order derivative 
can induce a spatial order, because the energy of the system 
should have a minimum at a non-zero momentum, and the 

simplest such example is given by E{p) p^ + p'^- In this 

regard, we mention that many pattern-forming models in vari- 
ous fields have higher-order derivatives, such as the convective 
instability yffl, the magnetic fluids 13 ill , and the generalized 

GL theory PZI- 

Unfortunately, however, Eq. ©is not integrable in gen- 
eral. Thus the next problem is how to construct an integrable 
model including higher-order derivatives. This is solved as 
follows. The integrable model generally has infinitely many 
conserved quantities /i, /2, . . . , and the model of their linear 
combination is also solvable by the ISM. Since the even- 
number /„ in the NLS hierarchy breaks a parity symmetry 
I 22n . the minimal model including higher-order derivatives is 
given by // = C3/3 -I- C5/5, that is the model ([T]l. 

The system is unstable if C5 < since the dispersion of the 
linearized operator becomes e ~ -k^. We can also confirm 
that if both C3 and C5 are positive, the ground state is always a 
trivial uniform state. Thus, the non-trivial physics arises when 
C3 < and C5 > 0. Henceforth we consider this case. 

Density-modulated ground state. — Let us begin the analy- 
sis of the model ([T]) in detail. We first determine the static 
ground state. Although the most general stationary solution 
to Eq. (|2ll is a quasi-periodic solution written by the Riemann 
theta function with genus g - 3 11241 l33n . here we assume 
that higher genus solutions are energetically unfavored, and 
we only consider the following two candidates; the FF state 
and the LO state 

qFF{x,p,p)^ ^/pe'P\ (4) 

<]Lo{x, p,m)— yfma sn{ax\m), a := s/p/Qim), (5) 

where Q(m) := 1 - -^^ with K{m) and E{m) being the com- 
plete elliptic integral of the first and second kind, respectively 
I 34II . p is the average of particle number density, and p and 
m are variational parameters to be chosen to minimize the 
energy. We can check that these states solve Eq. (O if we 
set p - cjip^ + 2p) + csip'^ + \2pp^ + 6p") for q^^ and 
p - Ci{m+\)a^ + cs{m^ +Am+\)a'^ for q\^o. Let £ff(p, p) and 
£lo(p, m) be the energies per particle for FF and LO states, 
and let p - Pg(p) and m - nig{p) be the values which mini- 
mize £ff and £lo under fixed p, respectively. We also define 
Pg{m) as an inverse function of nigip). The evaluation of them 
is summarized in the Supplemental Material 13511 . Figure [T| 
shows the plot of £ff(p, Pgip)) and £lo(p, '«g(p)) and corre- 
sponding periods. From Fig.[T] we can observe the two facts: 
(i) If the number density becomes small (p < j^), the LO 
state has a lower energy than the uniform state. (ii)The LO 




FIG. I . (Color online) (upper) The energies per particle for FF and 
LO states. We set cj = -1 and C5 = 1. For reference, we also show 
the energy for the uniform state q = yfp. (lower) The periods of FF 
and LO states, which are given by 2;r/p and 4K/a, respectively. 



state always has a lower energy than the FF state, except for 
the dilute limit p — > 0. Thus the density-modulated LO state 
is shown to be the lowest-energy state for small particle den- 
sities. As shown later, the LO state is linearly stable. 

AKNS form. — Equation (O enjoys the Ablowitz-Kaup- 
Newell-Segur (AKNS) or the zero-curvature representation 
llllll]: 



dJ=U(x,t,A)f, d,f^V(x,t,A)f, 



(6) 



where /I is a spectral parameter and / is a two-component vec- 
tor. The matrices U and V for Eq. (O are given by l23ll2411 



U ■ 



-lA q 
q* i/l 



y = -jUy">H-C3V^^*-HC5y<^\ (7) 



where 0"^ = 2;":o(-2/l)"--'-'M<^' and the explicit forms of 
M^-'' are given in Ref. 1371 ] with setting r - q*. We can check 
that the compatibility condition U, - Vx + [U,V] = yields 
Eq. (O. Note that the first equation of Eq. Q is equivalent 
to the linearized fermionic BdG equation for the quasiparticle 
wavefunctions with eigenenergy -A: 



-id, A 
A* 15, 



/ = -Af, A = iq. 



(8) 



Let us define V'3 :- V|c3=i, e5=0' which yields the AKNSi 
equation iq, = -pq - q,, + 2\q\^q. Using V'3, the associated 
Riemann surface for the LO state (|5]l is given by 

w' = det Vilg^,,^^ = 4/ - 20-2(1 + m)A^ + ^a^il - mf. (9) 

The condition or > determines the continuous spectrum 
of the BdG operator i37ll . which is given by \A\ < ^ 

and \A\ > 2 ^^^^ two-gap structure has been known in 
both the condensed matter and the algebro-geometric context 
lUmm. We also define 



o]:- 



det yL=„ ^ = [AcsA^ + C3 + csa\m + 1)]-^^ (10) 
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FIG. 2. Relation between the uniformization variable z and the spec- 
tral parameter -A. Here, -A can be interpreted as the eigenenergy 
of the fermionic BdG operator [Eq. {S}]. Regions with the same cir- 
cled numbers correspond to each other. The scattering states exist 
on the lines Imz = and \mz = K. The cross marks on the lines 
Rez = ±K'/2 represent discrete eigenvalues for bound states. In this 
example there are five bound states, so it represents a five-soliton so- 
lution. The rectangular contour shown by the bold line is used to 
derive the GLM equation l ll7t . 



which will become necessary to describe the time evolution of 
transition coefficients. 

Eigenstates of fermionic BdG operator. — We parametrize 
the Riemann surface (|9]) by the uniformization variable z'. 



-Kz) = 



\a cn(i2|m) dn(i2|m) 
2 sn(iz|m) 

„2 



(11) 



w(z) = aA'iz) = — [dn^acr - z)\m) - dn^CizIm)] . (12) 

The relation between z and -A is shown in Fig. [31 ws is 
also parametrized in the same way: u)s(z) - [4cg/l(z)^ -H C3 -)- 
c=,a^{m+ l)](xt(z). A{z) and w(z) have the following periodicity, 
parity, and the involution relation: A{z) - Mj, + InK' H- 7\IK) - 
A(K'-z) = AizT and w(z) = Cii{z+2nK'+2UK) = -a)(K'-z) = 
cl)(z*)* with / and n integers, a)s(z) has the same property with 

Henceforth we use the shorthand notation z' - K' -z unless 
otherwise noted. The two linearly independent solutions for 
Eq. ^ with A - A{z) and q - q-Lo{x,p, m) are given by 



faix, z) 



iQ;j?2(0)T?4(0)e' 



\ka(z)x 



??i(X-iZ)/i?4(iZ) 
)?3(0))?4(^) \^a(X - iZ)l-9i (iZ) 



(13) 



and f{){x, z'), where ■fiaiu) — '&a(u, q) is a theta function with 
the nome q = q-^^'Ik^ ^ ^ W'^ = #' and hiz) = k(z) - 



^ with a crystal momentum kiz) := ^(z)(l - M^n-f^M) + 

t(L^F^J -h |j. Here n{n\m) is the complete elliptic integral 
of the third kind and L»J is a floor function 13411 . Note that k{z) 
is meromorphic in the whole z-plane, since the discontinuity 
of n and L»J cancels out. k(z) satisfies k(z) = k{z + 2nK' + 
2UK) -H 2f = -k(z') = k(zy. /o(x,z) has the periodicity 
/o(x, z) = (- l)'/o(jc, z+2nA''+2i/^) and the involution relation 
fo(x,z') - o-ifo(x,z*y ■ Note that the solution ( fTSl l is also 




FIG. 3. (Color online) Bosonic Bogoliubov spectrum for the sn state 
[Eq. Q]. Here feose, <?B ose) = (2fc(z), -Icosiz))- We set C3 = -1, cs = 
\,m = 0.3, and a = -^Pg(in)/Q{m) = 0.638. The black solid line 
represents the branch of Bogoliubov phonons and corresponds to ® 
and (3) in Fig. [2] The red dashed line represents the branch of lattice- 
vibrating phonons and corresponds to (3) in Fig. [2] 



found in Ref. 14011 with a slightly different convention. 

The time evolution of /o(x, z) by Eq. ^ with A - A(z) is 
described using W5. If we choose the initial condition at f = 
as / = A/o(x, z) + Bfo{x,z'), the state at t is given by / = 
Afoix, z)e''^'fe)' -t- Bfo(x, z')e-'"'^'^'. 

Eigenstates of bosonic Bogoliubov operator and linear 
stability. — In the previous paragraph, we have obtained the 
spectrum of the fermionic BdG operator by regarding ^ as a 
fermionic condensate. Next, let us derive the bosonic Bogoli- 
ubov spectrum by regarding ^ as a bosonic condensate. The 
bosonic Bogoliubov equation can be obtained by linearization 
of E g. Q | 4in. and we can solve it by the squared eigenfunc 
tion 1421 14311 . More specifically, let / = (Mpermi, vpermi)^ 

the solution of Eq. ^, and (Mbosc, Vbosb) = ("Fermi' ~^Fenni) 

solves the bosonic Bogoliubov equation. See the Supplemen- 
tal Material for technical details 13511 . Using this fact and the 
explicit solution ( fT3b . we can readily obtain the dispersion re- 
lation as (2A:(z), -2w5(z)). See Fig. [3] Since this bosonic con- 
densate breaks the U(l)-gauge and the translational symmetry, 
two NG modes appear, i.e., the Bogoliubov phonon and the 
lattice- vibrating phonon. Figure [3] also automatically proves 
that there is no negative or complex eigenvalue. Thus the LO 
state is stable Il44tl . 

ISM with LO background. — From this point onwards, we 
formulate the ISM in the presence of the LO background, and 
construct the soliton solution. The corresponding problem for 
the KdV equation is solved in Ref. [45n and generalized to the 
finite-gap potential case in Ref. 1461) . Here we construct the 
corresponding theory for the NLS hierarchy. In this Letter we 
only show the main result and sketch of derivations. Detailed 
discussions will be published elsewhere. 

Let us assume that q satisfies the boundary condition 



be 



q{x) 



q-Loix) 

qLo(x - .x:o)e2i^o 



(x — > —00) 
(x — > -1-00), 



(14) 



where we omit the arguments p and m. (po and xq denote the 
phase and density shift due to solitons and ripple waves. We 
define the right Jost solution by the asymptotic form 



/+(x,z) 



fl(z)/o(x, z) + b{z)fa{x, z') {x ■ 



-00) 



-Vo(x-xo,z) 



(x — > -1-00). 



(15) 




FIG. 4. (Color online) Relation between the velocity of the soliton 
v and the spectral parameter -A when m ^ 1, where m is an ellip- 
tic parameter. The shaded areas represent the regions of continuous 
spectra. The soliton with zero velocity is a static dislocation. The 
other two represent the dark soliton billiard (Figure [Ha)) and the 
gray solitons passing through the dark-soliton array with a smaller 
interaction. Note that this figure is exaggerated and the continuous 
spectrum around /i = is very narrow when m ^ 1. 



This expression also defines the transition coeflicients a{z) and 
b{z). The other right Jost solution is given by /+(x, z')- Cal- 
culating the Wronskian W - f+{x,z'Y (}cr2)f+{x,z), the rela- 
tion a(z)a{z') - b{z)b{z') - 1 follows. The left Jost solution 
is defined by /_(x, z) — > /o(.x:, z')- The right and left Jost 
solutions are related as follows: 

(/,(x,z) A(x,z')) = (/-(x,z') f-i^^z))i^^% ^JJ5)- (16) 

f± has the same periodicity and involution relation with /q. 
The transition coefficients satisfy a{z') - a{z*)* and b{z') - 
b(z*y and have the same periodicity with A{z). 

Discrete eigenvalues for bound states are given by the ze- 
ros of a(z). Since the BdG operator is self-adjoint, bound 
states appear only for real A, which corresponds to Re z = ±i, 
as shown in Fig. [J] Let us assume that there are N zeros 
and write them as zi, . . . , Za^. We can prove that the normal- 
ization constant for the y'-th bound state defined by Cy^ := 

n dx/_(x, Zj^'f-ix, Zj) is given by Cf = 2iaa'izj)/bizj). 

We introduce the triangular representation for the left Jost 
solution: f-{x,z) - fo{x,z') + J_^dxr(x,y)fo(y, z')- Here the 
integral kernel r{x,y) is a 2 x 2 matrix. We can construct 
q(x) using this kernel: q(x) - qi^oix) + 2ri2(x, x). Calculat- 
ing J^ dz^[the first column of Eq. ^-foix, z)]fo(y, z'fcTi, 
where /? is a rectangular contour shown in Fig. |2l we obtain 
the Gel'fand-Levitan-Marchenko (GLM) equation 



r(x,y) + Q( 






AwT{x, w)^{w,y) -0, iy < x). 



Q(x,y):= I -—fQ{x,z')fQ{y,z'fcri 

Jr 4n a(z) 

+ aZ%iCJMx,z'j)My,z'/cTi (17) 

withz' :- K' - Zi- 

J ^ 

N-soliton solution. — Let us solve the GLM equation for the 
reflectionless case b{z) - 0. It can be solved by the following 
ansatz: T{x,y) = T!^^^Cj[l'l^)fo{y,z:p'(Ti. Substitution of 
this ansatz into the GLM equation yields the linear equations 




FIG. 5. (Color online) Examples of the one-soliton solution 
[Eq. M9V\. The plot shows the amplitude \q{x,i)\'-. (a) Dark soli- 
ton billiard. The parameters are c^ = -1, C5 = 1, m = 0.999, a = 
Vp«(™)/2('") = 0.527, zi = -0.5r + 03iK (<^ A = 0.0576), and 
Ci(0) = 0.160. The velocity is vi = 0.0711. (b) Snap shot of the 
envelope soliton when the background sn function is almost trigono- 
metric. The arrow shows the direction of the soliton propagation. 
The parameters are C3 = -1, C5 = 1, ra = 0.3, a = 0.638, zi = 
-0.5K' + 0.55iii: {■^ A = 0.243), Ci(0) = 5.49, and t = 40. The 
velocity is vi = 0.239. See also the ancillary animated GIF files. 



for gjixYs: 

gj(x) + aCjUo(x, z'j) + aZli CjCiM(x, Zj, zi)giix) = 0, (18) 
where uq(x, z) is the first component of /o(x, z) and 



Mix,Zj,Zi) := -20- 






i?3(0)7?4(^)»?i(iZ;-iz;) 



with Zj - -j^ and Z^' = 



"^^j^"'^ ■ Using g/(A:)'s, we obtain 
the A'-soliton solution q{x) - qi,o{x) + lYijCjgj{x)UQ{x,z!^. 
The time-dependence of C; is given by Cj(t) = Cj(0)e"''^=(^')' 



and the velocity of each soliton becomes v j - - ^^"^, ^^, . 
particular, the one-soliton solution is given by 



In 



9(x, = q\^o{x) ■ 



laC^iW^ 



2p-2ia)5(C])/ 



uq{x,K' -zi) 



1 -H Q'Ci(0)2e-2'"5bi)'M(x,zi,Zi) '' 



(19) 



This one-soliton solution shows various behaviors dependent 
on the choice of parameters. When m ^ 1, we can broadly 
classify it into three categories by its velocity, that is, dark 
soliton billiards, stationary dislocations, and gray solitons. 
See Fig. |4l In this case, the soliton propagation can be un- 
derstood approximately as a successive collision between the 
moving soliton and the array of static dark solitons. Fig- 
ure 13 a) shows an example of the dark soliton billiard. We 
note that the LO background experiences a density shift af- 
ter the passing of the soliton, as opposed to conventional gap 



solitons. When m ^ 0, the LO background is almost trigono- 
metric. In this case the distinction between bilhards and gray 
solitons becomes obscure, and solitons with any velocity can 
be observed as envelope solitons (Figure[5jb)). This behavior 
is very similar to the solitons observed in the soft-core bosons 
12911 ■ See also the ancillary animated GIF files for the other ex- 
amples. The soliton parameters for these animations are noted 
in the Supplemental Material 1351]. 

Summary. — In summary, we have introduced an integrable 
model of density-modulated quantum condensates and have 
provided an A'-soliton solution. We have shown that the model 
can be constructed as a linear combination of conserved quan- 
tities in the NLS hierarchy and that density-modulated states 
have lower energy than uniform states. The obtained exact 
soliton solutions contain various novel dynamics, such as soli- 
ton billiards, stationary dislocations, gray solitons, and enve- 
lope solitons. We believe that our results are universal and 
will be useful to understand nonequilibrium and transport phe- 
nomena in non-uniform quantum matters. Finally we note 
several future works, such as the extension to the self- focusing 
case (r - -q*), the determination of self-consistent soliton 



solutions in BdG systems 114711 . and an exact theory of spon- 
taneous translational-symmetry breaking in the corresponding 
quantum model. 
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Supplemental Material 



I. ENERGY EVALUATION FOR FF AND LO STATES 

The energy density h(x) at a position x is defined by the in- 
tegrand of Eq. ^ of the main article. The energy per particle 
is then defined hy & - J dxli(x) J dx\if/\^, where L is a pe- 
riod and given by L = 2;r//7 for the FF state and L = 4K/a for 
the LO state, respectively. After a little calculation, we obtain 



£ff(p, p) = C3(p^ +p) + csip'^ + 6p-p + 2p-), 
p[m + (m + l)Q{m)] 



Slo(p, m) = C3 



+ C5- 



3G(m)2 
p^[2m(m +1) + {jn^ + 4m + l)Q(m)] 



(SI) 



(S2) 



Let p - Pg(p) and m - ing(p) be the value which minimize the 
above energy density fipF and £lo with fixed p, respectively. 
Then, they are determined as follows: 



Pg(p) 



nigip) 







/I 



: V-(C3 + 6C5P)/(2C5) {p<^), 

(p>t|) 

[inverse function of pg(m) (p < TgT')' 

-5c3[-2m + {\+ m)Q(m)]Q(m) 
6c5[-3m(l +m) + (l +4m + m^)Q(m)] ' 

Here we have assumed C3 < and C5 > 0. 



Pg(m) := 



(S3) 

(S4) 
(S5) 



II. SOLUTIONS FOR BOSONIC BOGOLIUBOV EQUATION 

The bosonic Bogoliubov equation can be obtained by lin- 
earizing Eq. (|2]) of the main article. More specifically, setting 
q = q + 6q in Eq. Q and cutting higher order terms with 
respect to Sq, we obtain 

idq, = -pSq + c^l-dq,,;, + 2{2\q\^5q + q^dq*)] + C5[6q,,xxx 

- 2{q*6q + q6q*)xxq - 2{\qf)xxSq - 3(q^)xxSq* 

- 6q*iq6q)xx + 6\q\H3\q\^6q + 2q^6q*)]. (S6) 

and c.c. of 



Substituting {6q, 6q*) - (m, -v) into Eq. 
Eq. (IS6b . we obtain the Bogoliubov equation: 

m, - -fiu + C3[-Uxx + 2(2\q\^u - q^v)] + Csluxxx 

- 2(q*u - qv)xxq - 2(\qf)xxU + 3(^^).„v 

- 6q*(qu)xx + 6\q\\3\qfu - 2q^v)l 

h'l = ^v - C3[-v.„ + 2(2\qfv - q*^u)] - Ci[vxxxx 

- 2(qv - q*u)xxq* - 2(l9l').„v + 3(q*^)xxU 
-6qiq'v)xx + 6\q\\3\q\'v-2q*^u)]. 



(S7) 



(S8) 



The stationary Bogoliubov equation with the eigenenergy 
e can be obtained by the substitution {u(x, t), v(x, t)) - 
(m(x), v(x))e""^'. As stated in the main article, we can solve 
the above equation by the squared eigenfunction. If / = 



(Mpermi, vpermi)^ is a solution of Eq. Q of the main article. 



("Bose- l^Bose) - («; 



2 
Fermi ' 



Fermi 



(S9) 



becomes a solution of the bosonic Bogoliubov equation 
dSTI l and JSSI l. Since the time-evolution of the solution 
("Fermi, Vpermi) wlth the spectral parameter A - A{z) is given by 
giw5(z)'^ the time-dependence of the corresponding bosonic so- 
lution is given by e^™'*'*', which means that the eigenenergy 
of the eigenstate labeled by z is given by - 2^5(2). For the 
same reason, the crystal momentum is given by 2k{z). Thus 
we obtain Fig. [3] of the main article. 



III. PARAMETERS FOR THE ANIMATED GIF FILES 

• animationl.gif: Dark soliton billiard. 
The parameters are the same as Fig.jS^a). 

• animation2.gif: Gray soliton. 

zi = -0.5K' + 0.05iK (« ^ = 0.471) and Ci(0) = 
1.23. The velocity is given by vi = -0.469. The other 
parameters are the same as Fig. [3 a). 

• animation3.gif: Static dislocation. 

zi = -Q.5K' + 0.\Q66\K (^ ^ = 0.333) and Ci(0) = 
0.0463. The velocity is given by vi — 0. The other 
parameters are the same as Fig. 13 a). 

• animation4.gif: Example of 3-soliton solution. 

zi = -Q.5K' + 0.3iK, zi = -0.5^' + 0.l066iK, Z3 = 
-Q.5K' + 0.051/*:, Ci(0) = 0.643, C2(0) = 0.0463, and 
C3(0) = 243. The other parameters are the same as 
Fig. Ha). 

• animation5.gif: Envelope soliton. 

The parameters are the same as Fig.|3b). 

• animation6.gif: Another envelope soliton. 

zi = -0.5K' + 0.341^ (« i = 0.357) and Ci(0) = 
5.17. The velocity becomes vi = -0.0477. The other 
parameters are the same as Fig.|5jb). 



